Practical 1B - microbiota data

NUTRIM microbiome & metabolome workshop

Author

David Barnett

Published

May 27, 2024

1 Intro

Now we’re going to look at the microbiota data.

These data have already been processed into a table, counts of how often each sequence occurs in each sample.

It started as huge fastq text files, full of As, Cs, Ts and Gs, but we will not practice sequence data processing today, because it takes quite a long time to run.

If you generate your data with us at Medical Microbiology, we will run the latest processing approaches on our computational infrastructure and provide you with the processed data.

Nowadays we do amplicon sequencing with Illumina MiSeq, HiSeq, or similar technologies, and denoise the output with DADA2 to produce ASV count tables.

The example data we will use are a little older. The amplicons were sequenced with 454 pyrosequencing and clustered into OTUs, but the core principles of analysis remain the same.

And as you will see later today, the same approaches can also be applied to taxonomic abundance count tables obtained from shotgun metagenomic sequencing.

First we’re going to use the R skills we practised in part A to inspect this data.

  • In this practical, solution blocks like this are available…

  • But try and write code independently before looking at the answer!

  • Attempting to recall knowledge or solve problems is proven to enhance learning.

  • So don’t look unless you’re really stuck! 💪

After that, we’ll take a look at some specialist R packages for microbiome analysis.

2 Learning goals 💪

3 Load R packages 📦

here() starts at /Users/david/Documents/git/R-projects/teaching/workshops/2024-NUTRIM-microbiome
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

4 Read and inspect data 🔍

Read the metadata file from part 1A, using the RDS version.

meta <- read_rds(file = here("data/papa2012/processed/all_metadata.rds"))

Read the count table: this is stored as a TSV (tab-separated variables) formatted text file.

counts <- read_tsv(file = here("data/papa2012/papa2012_OTU_count_table.txt"))
Rows: 91 Columns: 36350
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr     (1): sample
dbl (36349): OTU_00001, OTU_00002, OTU_00003, OTU_00004, OTU_00005, OTU_0000...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
counts
# A tibble: 91 × 36,350
   sample OTU_00001 OTU_00002 OTU_00003 OTU_00004 OTU_00005 OTU_00006 OTU_00007
   <chr>      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 099_AX         0        56        30         0         0       103         0
 2 199_AX        29         0         9        17         4         5        57
 3 062_BZ         0        95        33         0         0         0        32
 4 194_AZ        75        49       228         0        30        35         0
 5 166_AX       138         3        12         0        13        88       103
 6 219_AX         0       347         2         2         0         0         1
 7 132_AX       121         0         7         0         0         0       181
 8 026_AX         0         0         0         0         0         0         0
 9 102_AZ         0         0         0         0         0         0         0
10 140_AX        47         0        16         1         1        46         4
# ℹ 81 more rows
# ℹ 36,342 more variables: OTU_00008 <dbl>, OTU_00009 <dbl>, OTU_00010 <dbl>,
#   OTU_00011 <dbl>, OTU_00012 <dbl>, OTU_00013 <dbl>, OTU_00014 <dbl>,
#   OTU_00015 <dbl>, OTU_00016 <dbl>, OTU_00017 <dbl>, OTU_00018 <dbl>,
#   OTU_00019 <dbl>, OTU_00020 <dbl>, OTU_00021 <dbl>, OTU_00022 <dbl>,
#   OTU_00023 <dbl>, OTU_00024 <dbl>, OTU_00025 <dbl>, OTU_00026 <dbl>,
#   OTU_00027 <dbl>, OTU_00028 <dbl>, OTU_00029 <dbl>, OTU_00030 <dbl>, …

Your first challenge: read the taxonomy table stored in "data/papa2012/papa2012_taxonomy_table.txt" (call the object taxonomy)

taxonomy <- read_tsv(file = here("data/papa2012/papa2012_taxonomy_table.txt"))
Rows: 36349 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): OTU, Kingdom, Phylum, Class, Order, Family, Genus

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
taxonomy
# A tibble: 36,349 × 7
   OTU       Kingdom  Phylum         Class               Order      Family Genus
   <chr>     <chr>    <chr>          <chr>               <chr>      <chr>  <chr>
 1 OTU_00001 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 2 OTU_00002 Bacteria Proteobacteria Gammaproteobacteria Enterobac… Enter… Esch…
 3 OTU_00003 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Blau…
 4 OTU_00004 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Rose…
 5 OTU_00005 Bacteria Bacteroidetes  Bacteroidia         Bacteroid… Bacte… Bact…
 6 OTU_00006 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 7 OTU_00007 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Rose…
 8 OTU_00008 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 9 OTU_00009 Bacteria Proteobacteria Gammaproteobacteria Enterobac… Enter… Esch…
10 OTU_00010 Bacteria Firmicutes     Negativicutes       Selenomon… Veill… Dial…
# ℹ 36,339 more rows

4.1 Taxonomy table

Now practice inspecting the taxonomy table by completing the following tasks:

Check how many distinct genera there are. Tip: use unique() and length()

taxonomy$Genus %>%
  unique() %>%
  length()
[1] 142

How many OTUs are there in each Phylum? Tip: use count() or table()

taxonomy %>% count(Phylum)
# A tibble: 9 × 2
  Phylum                        n
  <chr>                     <int>
1 Actinobacteria              307
2 Bacteroidetes             13109
3 Cyanobacteria/Chloroplast     1
4 Firmicutes                17653
5 Fusobacteria                275
6 Proteobacteria             4806
7 Synergistetes                 6
8 Verrucomicrobia              59
9 <NA>                        133
taxonomy$Phylum %>% table(useNA = "ifany")
.
           Actinobacteria             Bacteroidetes Cyanobacteria/Chloroplast 
                      307                     13109                         1 
               Firmicutes              Fusobacteria            Proteobacteria 
                    17653                       275                      4806 
            Synergistetes           Verrucomicrobia                      <NA> 
                        6                        59                       133 

What genera are in the phylum Actinobacteria? Tip: use filter()

taxonomy %>%
  filter(Phylum == "Actinobacteria") %>%
  count(Genus, sort = TRUE)
# A tibble: 9 × 2
  Genus               n
  <chr>           <int>
1 Bifidobacterium   170
2 Collinsella       111
3 Eggerthella         8
4 Actinomyces         7
5 <NA>                5
6 Gordonibacter       3
7 Asaccharobacter     1
8 Atopobium           1
9 Rothia              1

Or for just their names:

taxonomy %>%
  filter(Phylum == "Actinobacteria") %>%
  pull(Genus) %>%
  unique()
[1] "Collinsella"     "Bifidobacterium" "Eggerthella"     "Atopobium"      
[5] "Actinomyces"     NA                "Gordonibacter"   "Asaccharobacter"
[9] "Rothia"         

4.2 OTU count table

First let’s plot a histogram of OTU number 1.

counts %>% ggplot(aes(OTU_00001)) +
  geom_histogram(binwidth = 5)

Looks like there are a lot of zeros!

table(OTU00001_has_0_counts = counts$OTU_00001 == 0, useNA = "ifany")
OTU00001_has_0_counts
FALSE  TRUE 
   59    32 

Attempt the following tasks to explore further!

Use filter() to plot only the non-zero entries for OTU 1.

Try also using + scale_x_log10() to transform the plot axis scale.

counts %>%
  filter(OTU_00001 != 0) %>%
  ggplot(aes(OTU_00001)) +
  geom_histogram(bins = 30) +
  scale_x_log10()

Use mutate() to create a temporary log-transformed variable log_OTU1 and plot its distribution.

Note that you can’t do log(0), so you will need to add 1 to all values first, i.e. log10(OTU_00001 + 1)

counts %>%
  mutate(log_OTU1 = log10(OTU_00001 + 1)) %>%
  ggplot(aes(log_OTU1)) +
  geom_histogram(bins = 30)

Plot OTU 1 against OTU 2 as a scatter plot using geom_point.

Remember to add 1 before log10 transformation of both variables.

counts %>%
  mutate(log_OTU1 = log10(OTU_00001 + 1)) %>%
  mutate(log_OTU2 = log10(OTU_00002 + 1)) %>%
  ggplot(aes(x = log_OTU1, y = log_OTU2)) +
  geom_point()

But what type of bacteria do these amplicon sequences belong to?

Look up the classifications of OTU 1 and OTU 2 in the taxonomy table.

taxonomy %>% filter(OTU %in% c("OTU_00001", "OTU_00002"))
# A tibble: 2 × 7
  OTU       Kingdom  Phylum         Class               Order       Family Genus
  <chr>     <chr>    <chr>          <chr>               <chr>       <chr>  <chr>
1 OTU_00001 Bacteria Firmicutes     Clostridia          Clostridia… Rumin… Faec…
2 OTU_00002 Bacteria Proteobacteria Gammaproteobacteria Enterobact… Enter… Esch…
taxonomy %>%
  filter(OTU %in% c("OTU_00001", "OTU_00002")) %>%
  select(OTU, Phylum, Family, Genus) 
# A tibble: 2 × 4
  OTU       Phylum         Family             Genus               
  <chr>     <chr>          <chr>              <chr>               
1 OTU_00001 Firmicutes     Ruminococcaceae    Faecalibacterium    
2 OTU_00002 Proteobacteria Enterobacteriaceae Escherichia/Shigella

Now make histograms and look up the taxonomy for the next thousand OTUs…

Okay, that was a joke. 🤡

It is possible to make a thousand plots, because R is very good at repetitive tasks.

But, this would not be very useful, because you could not look at them all.

In the next section of this practical, we will explore smarter ways to analyse microbiota compositions.

5 Assemble a phyloseq 🛠️

So far we have not used any R packages specialised for microbiome data.

Let’s do so, because it will make our next tasks a lot easier!

microViz version 0.12.1 - Copyright (C) 2021-2024 David Barnett
! Website: https://david-barnett.github.io/microViz
✔ Useful?  For citation details, run: `citation("microViz")`
✖ Silence? `suppressPackageStartupMessages(library(microViz))`

We will start by combining our three dataframes into one phyloseq object

5.1 OTU counts + taxonomy

phyloseq uses row or column names to match OTUs across the count and taxonomy tables.

In addition, the taxonomy table must contain only taxonomic ranks in descending order, and the count table must be converted to a numeric matrix.

count_matrix <- counts %>% select(-sample) %>% as.matrix()
rownames(count_matrix) <- counts$sample
tax_matrix <- taxonomy %>% select(!OTU) %>% as.matrix()
rownames(tax_matrix) <- taxonomy$OTU
ps <- phyloseq(
  otu_table(count_matrix, taxa_are_rows = FALSE), 
  tax_table(tax_matrix)
)
ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 91 samples ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

5.2 Adding sample metadata

phyloseq uses row names to match samples (across the OTU count table and the sample metadata).

[1] "099_AX" "199_AX" "062_BZ" "194_AZ" "166_AX" "219_AX"

Notice that in the phyloseq object the sample names have underscores, but they have hyphens in the sample metadata.

meta
# A tibble: 90 × 17
   ID    sample case_control diagnosis activity active   gender ethnicity
   <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr>  <chr>    
 1 099A  099-AX Case         UC        severe   active   female White    
 2 199A  199-AX Control      Other     control  control  female Other    
 3 062B  062-BZ Case         CD        mild     active   male   White    
 4 194A  194-AZ Case         UC        mild     active   male   White    
 5 166A  166-AX Case         UC        severe   active   female Black    
 6 219A  219-AX Case         UC        inactive inactive female <NA>     
 7 132A  132-AX Case         CD        mild     active   female White    
 8 026A  026-AX Case         UC        mild     active   male   White    
 9 102A  102-AZ Control      Other     control  control  male   White    
10 140A  140-AX Control      Other     control  control  female White    
# ℹ 80 more rows
# ℹ 9 more variables: family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>

We can fix that with the stringr package str_replace function.

meta$sample <- meta$sample %>% str_replace(pattern = "-", replacement = "_")
meta
# A tibble: 90 × 17
   ID    sample case_control diagnosis activity active   gender ethnicity
   <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr>  <chr>    
 1 099A  099_AX Case         UC        severe   active   female White    
 2 199A  199_AX Control      Other     control  control  female Other    
 3 062B  062_BZ Case         CD        mild     active   male   White    
 4 194A  194_AZ Case         UC        mild     active   male   White    
 5 166A  166_AX Case         UC        severe   active   female Black    
 6 219A  219_AX Case         UC        inactive inactive female <NA>     
 7 132A  132_AX Case         CD        mild     active   female White    
 8 026A  026_AX Case         UC        mild     active   male   White    
 9 102A  102_AZ Control      Other     control  control  male   White    
10 140A  140_AX Control      Other     control  control  female White    
# ℹ 80 more rows
# ℹ 9 more variables: family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>

Now make them row names, check they all match, and add the metadata to the phyloseq object.

meta_df <- as.data.frame(meta)
rownames(meta_df) <- meta$sample
all(rownames(meta_df) %in% sample_names(ps))
[1] TRUE
sample_data(ps) <- meta_df
ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 90 samples ]
sample_data() Sample Data:       [ 90 samples by 17 sample variables ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

6 microViz 🦠👁️

microViz provides tools for working with phyloseq objects.

Let’s take a look with some basic microViz functions.

samdat_tbl(ps) # retrieve sample data as a tibble
# A tibble: 90 × 18
   .sample_name ID    sample case_control diagnosis activity active   gender
   <chr>        <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr> 
 1 099_AX       099A  099_AX Case         UC        severe   active   female
 2 199_AX       199A  199_AX Control      Other     control  control  female
 3 062_BZ       062B  062_BZ Case         CD        mild     active   male  
 4 194_AZ       194A  194_AZ Case         UC        mild     active   male  
 5 166_AX       166A  166_AX Case         UC        severe   active   female
 6 219_AX       219A  219_AX Case         UC        inactive inactive female
 7 132_AX       132A  132_AX Case         CD        mild     active   female
 8 026_AX       026A  026_AX Case         UC        mild     active   male  
 9 102_AZ       102A  102_AZ Control      Other     control  control  male  
10 140_AX       140A  140_AX Control      Other     control  control  female
# ℹ 80 more rows
# ℹ 10 more variables: ethnicity <chr>, family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>
# get the OTU table, or part of it
otu_get(ps, taxa = 1:5, samples = c("132_AX", "166_AX", "102_AZ"))
OTU Table:          [5 taxa and 3 samples]
                     taxa are columns
       OTU_00001 OTU_00002 OTU_00003 OTU_00004 OTU_00005
132_AX       121         0         7         0         0
166_AX       138         3        12         0        13
102_AZ         0         0         0         0         0
# get the taxonomy table
tt_get(ps) %>% head(3)
Taxonomy Table:     [3 taxa by 6 taxonomic ranks]:
          Kingdom    Phylum           Class                 Order              
OTU_00001 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales"    
OTU_00002 "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Enterobacteriales"
OTU_00003 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales"    
          Family               Genus                 
OTU_00001 "Ruminococcaceae"    "Faecalibacterium"    
OTU_00002 "Enterobacteriaceae" "Escherichia/Shigella"
OTU_00003 "Lachnospiraceae"    "Blautia"             

7 Build some bar charts 📊

Sequencing data are compositional. The total number of reads per sample is mostly arbitrary, and so the counts should be interpreted as relative abundances instead of absolute abundances.

A simple way to visualise compositional data is as percentages, proportions of a whole.

Stacked bar charts are great for this: each bar represents one sample, and we can show the abundance data as proportions, after aggregating the counts by taxonomy.

ps %>% tax_fix() %>% comp_barplot("Phylum", n_taxa = 4, label = NULL)

This plot is aggregated into phyla, which is easy to read, and provides basic information.

We see that most samples contain a mixture of Firmicutes and Bacteroidetes, but some appear to be dominated by Proteobacteria instead.

Beware, phylum names have changed recently!

  • Actinobacteria is now Actinomycetota
  • Bacteroidetes is now Bacteroidota
  • Proteobacteria is now Pseudomonadota (!)
  • Firmicutes is now Bacillota (!!)

In all prior research, you will see the original names, but in coming years, the new names will be increasingly adopted.

The best source for checking/searching official and alternative names is probably LPSN at bacterio.net

ps %>% tax_fix() %>% comp_barplot("Genus", n_taxa = 11, merge_other = F, label = NULL)

This plot is aggregated into genera, which provides more detail, but is harder to read, and we cannot give every genus a distinct colour.

We see that many samples contain a relatively large proportion of Bacteroides or Faecalibacterium but there is quite some further variation!

8 Discover diversity 🌳🌲🌴

As a last task for this introductory session, we will calculate and visualise alpha diversity.

We will calculate Shannon diversity at the level of Genus.

  • The Shannon index is a commonly used diversity measure, with this formula: \(H = -\sum_{i=1}^Np_i\ln p_i\)
  • Shannon index is often is denoted with \(H\), and here \(p\) denotes proportion of the \(i\)’th taxon.
  • For each taxon \(i\), you multiply its proportional abundance \(p_i\) by the natural log of that proportion \(\ln p_i\), and sum those values.
  • Try it out for yourself to convince yourself you get larger (negative) values for higher proportions.
  • The highest value you can achieve with \(N\) taxa occurs with equal proportions
  • e.g. with 20 taxa, maximum diversity occurs if each has a proportion of 5%
  • Lastly, you change the sign of the result to a positive number, for ease of interpretation
  • This just makes more intuitive sense: as higher positive numbers indicates higher diversity.

We will also compute a exponentiated version, the effective number of genera, or effective Shannon index.

The exponent of the Shannon index \(e^H\) represents the number of taxa (genera) that would be present in an evenly abundant ecosystem with the same Shannon index.

  • The numeric value of the Shannon index itself has no intuitive meaning.
  • You can compare them, but can’t easily interpret any one number.
  • So, the concept of “effective numbers” of taxa is useful here.
  • If your original ecosystem was actually perfectly even, then \(e^H = N\)
  • Where N is the observed richness
  • The more uneven an ecosystem, the further \(e^H\) will be from \(N\)
ps_alpha <- ps %>% 
  tax_fix() %>%
  ps_calc_diversity(index = "shannon", rank = "Genus", varname = "Shannon_Genus") %>% 
  ps_mutate(Effective_Shannon_Genus = exp(Shannon_Genus))

ps_alpha
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 90 samples ]
sample_data() Sample Data:       [ 90 samples by 19 sample variables ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

The new genus-level Shannon diversity variables are stored in the sample data slot of ps_alpha

ps_alpha %>% samdat_tbl() %>% select(sample, Shannon_Genus, Effective_Shannon_Genus)
# A tibble: 90 × 3
   sample Shannon_Genus Effective_Shannon_Genus
   <chr>          <dbl>                   <dbl>
 1 099_AX          1.89                    6.63
 2 199_AX          2.71                   15.1 
 3 062_BZ          1.66                    5.25
 4 194_AZ          1.52                    4.59
 5 166_AX          1.92                    6.81
 6 219_AX          1.22                    3.39
 7 132_AX          2.35                   10.5 
 8 026_AX          2.14                    8.47
 9 102_AZ          2.56                   12.9 
10 140_AX          1.32                    3.74
# ℹ 80 more rows
ps_alpha %>% samdat_tbl() %>% pull(Shannon_Genus) %>% hist(main = "Shannon Diversity, Genus")

ps_alpha %>% samdat_tbl() %>% pull(Effective_Shannon_Genus) %>% hist(main = "Effective no. genera")

8.1 Interpreting diversity

Now each sample is labelled with the effective Shannon diversity of genera (\(e^H\)) - do you see a relationship to sample composition?

Code
ps_alpha %>% 
  ps_arrange(Effective_Shannon_Genus) %>% 
  ps_mutate(short_shannon = formatC(Effective_Shannon_Genus, digits = 2, format = "f")) %>% 
  comp_barplot(
    tax_level = "Genus", n_taxa = 19, merge_other = FALSE, 
    sample_order = "asis", label = "short_shannon"
  ) +
  coord_flip()

9 Next! ⏩

  • Its time to stop exploring and take a break.
  • In the next lecture you will learn more about the main approaches for microbiota data analysis.
  • In the next practical session, we will come back to these data with a plan, an analysis plan!

10 Session info

sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.0 (2024-04-24)
 os       macOS Sonoma 14.5
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Amsterdam
 date     2024-05-27
 pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version  date (UTC) lib source
 ade4               1.7-22   2023-02-06 [1] CRAN (R 4.4.0)
 ape                5.8      2024-04-11 [1] CRAN (R 4.4.0)
 Biobase            2.64.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 BiocGenerics       0.50.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 biomformat         1.32.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 Biostrings         2.72.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 bit                4.0.5    2022-11-15 [1] CRAN (R 4.4.0)
 bit64              4.0.5    2020-08-30 [1] CRAN (R 4.4.0)
 ca                 0.71.1   2020-01-24 [1] CRAN (R 4.4.0)
 cellranger         1.1.0    2016-07-27 [1] CRAN (R 4.4.0)
 cli                3.6.2    2023-12-11 [1] CRAN (R 4.4.0)
 cluster            2.1.6    2023-12-01 [1] CRAN (R 4.4.0)
 codetools          0.2-20   2024-03-31 [1] CRAN (R 4.4.0)
 colorspace         2.1-0    2023-01-23 [1] CRAN (R 4.4.0)
 crayon             1.5.2    2022-09-29 [1] CRAN (R 4.4.0)
 data.table         1.15.4   2024-03-30 [1] CRAN (R 4.4.0)
 digest             0.6.35   2024-03-11 [1] CRAN (R 4.4.0)
 dplyr            * 1.1.4    2023-11-17 [1] CRAN (R 4.4.0)
 evaluate           0.23     2023-11-01 [1] CRAN (R 4.4.0)
 fansi              1.0.6    2023-12-08 [1] CRAN (R 4.4.0)
 farver             2.1.2    2024-05-13 [1] CRAN (R 4.4.0)
 fastmap            1.2.0    2024-05-15 [1] CRAN (R 4.4.0)
 forcats          * 1.0.0    2023-01-29 [1] CRAN (R 4.4.0)
 foreach            1.5.2    2022-02-02 [1] CRAN (R 4.4.0)
 generics           0.1.3    2022-07-05 [1] CRAN (R 4.4.0)
 GenomeInfoDb       1.40.1   2024-05-24 [1] Bioconductor 3.19 (R 4.4.0)
 GenomeInfoDbData   1.2.12   2024-05-26 [1] Bioconductor
 ggplot2          * 3.5.1    2024-04-23 [1] CRAN (R 4.4.0)
 glue               1.7.0    2024-01-09 [1] CRAN (R 4.4.0)
 gtable             0.3.5    2024-04-22 [1] CRAN (R 4.4.0)
 here             * 1.0.1    2020-12-13 [1] CRAN (R 4.4.0)
 hms                1.1.3    2023-03-21 [1] CRAN (R 4.4.0)
 htmltools          0.5.8.1  2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets        1.6.4    2023-12-06 [1] CRAN (R 4.4.0)
 httr               1.4.7    2023-08-15 [1] CRAN (R 4.4.0)
 igraph             2.0.3    2024-03-13 [1] CRAN (R 4.4.0)
 IRanges            2.38.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 iterators          1.0.14   2022-02-05 [1] CRAN (R 4.4.0)
 jsonlite           1.8.8    2023-12-04 [1] CRAN (R 4.4.0)
 knitr              1.46     2024-04-06 [1] CRAN (R 4.4.0)
 labeling           0.4.3    2023-08-29 [1] CRAN (R 4.4.0)
 lattice            0.22-6   2024-03-20 [1] CRAN (R 4.4.0)
 lifecycle          1.0.4    2023-11-07 [1] CRAN (R 4.4.0)
 lubridate        * 1.9.3    2023-09-27 [1] CRAN (R 4.4.0)
 magrittr           2.0.3    2022-03-30 [1] CRAN (R 4.4.0)
 MASS               7.3-60.2 2024-04-24 [1] local
 Matrix             1.7-0    2024-03-22 [1] CRAN (R 4.4.0)
 mgcv               1.9-1    2023-12-21 [1] CRAN (R 4.4.0)
 microbiome         1.26.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 microViz         * 0.12.1   2024-05-03 [1] https://david-barnett.r-universe.dev (R 4.4.0)
 multtest           2.60.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 munsell            0.5.1    2024-04-01 [1] CRAN (R 4.4.0)
 nlme               3.1-164  2023-11-27 [1] CRAN (R 4.4.0)
 permute            0.9-7    2022-01-27 [1] CRAN (R 4.4.0)
 phyloseq         * 1.48.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 pillar             1.9.0    2023-03-22 [1] CRAN (R 4.4.0)
 pkgconfig          2.0.3    2019-09-22 [1] CRAN (R 4.4.0)
 plyr               1.8.9    2023-10-02 [1] CRAN (R 4.4.0)
 purrr            * 1.0.2    2023-08-10 [1] CRAN (R 4.4.0)
 R6                 2.5.1    2021-08-19 [1] CRAN (R 4.4.0)
 Rcpp               1.0.12   2024-01-09 [1] CRAN (R 4.4.0)
 readr            * 2.1.5    2024-01-10 [1] CRAN (R 4.4.0)
 readxl           * 1.4.3    2023-07-06 [1] CRAN (R 4.4.0)
 registry           0.5-1    2019-03-05 [1] CRAN (R 4.4.0)
 reshape2           1.4.4    2020-04-09 [1] CRAN (R 4.4.0)
 rhdf5              2.48.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 rhdf5filters       1.16.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 Rhdf5lib           1.26.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 rlang              1.1.3    2024-01-10 [1] CRAN (R 4.4.0)
 rmarkdown          2.27     2024-05-17 [1] CRAN (R 4.4.0)
 rprojroot          2.0.4    2023-11-05 [1] CRAN (R 4.4.0)
 rstudioapi         0.16.0   2024-03-24 [1] CRAN (R 4.4.0)
 Rtsne              0.17     2023-12-07 [1] CRAN (R 4.4.0)
 S4Vectors          0.42.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 scales             1.3.0    2023-11-28 [1] CRAN (R 4.4.0)
 seriation          1.5.5    2024-04-17 [1] CRAN (R 4.4.0)
 sessioninfo        1.2.2    2021-12-06 [1] CRAN (R 4.4.0)
 stringi            1.8.4    2024-05-06 [1] CRAN (R 4.4.0)
 stringr          * 1.5.1    2023-11-14 [1] CRAN (R 4.4.0)
 survival           3.5-8    2024-02-14 [1] CRAN (R 4.4.0)
 tibble           * 3.2.1    2023-03-20 [1] CRAN (R 4.4.0)
 tidyr            * 1.3.1    2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect         1.2.1    2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse        * 2.0.0    2023-02-22 [1] CRAN (R 4.4.0)
 timechange         0.3.0    2024-01-18 [1] CRAN (R 4.4.0)
 TSP                1.2-4    2023-04-04 [1] CRAN (R 4.4.0)
 tzdb               0.4.0    2023-05-12 [1] CRAN (R 4.4.0)
 UCSC.utils         1.0.0    2024-05-06 [1] Bioconductor 3.19 (R 4.4.0)
 utf8               1.2.4    2023-10-22 [1] CRAN (R 4.4.0)
 vctrs              0.6.5    2023-12-01 [1] CRAN (R 4.4.0)
 vegan              2.6-6.1  2024-05-21 [1] CRAN (R 4.4.0)
 vroom              1.6.5    2023-12-05 [1] CRAN (R 4.4.0)
 withr              3.0.0    2024-01-16 [1] CRAN (R 4.4.0)
 xfun               0.44     2024-05-15 [1] CRAN (R 4.4.0)
 XVector            0.44.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 yaml               2.3.8    2023-12-11 [1] CRAN (R 4.4.0)
 zlibbioc           1.50.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)

 [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library

──────────────────────────────────────────────────────────────────────────────